###############################################

## Weathering the Storm: Discordant Learning About Reputations for Reliability
## Replication for Scenario Based Experiment Results Article and Appendices

## Author: Bailee Donahue
## Last Updated: August 24, 2022


#################################################


rm(list=ls())
## Going to have to change the working directory
setwd("~/Downloads/dataverse_files")


library(MASS)
library(reshape2)
library(reshape)
library(countrycode)
library(states)
library(pltesim)
library(readstata13)
library(plyr)
library(foreign)
library(ggplot2)
##install.packages("cjoint")
library(cjoint)

##install.packages("FindIt")
library(FindIt)
#install.packages("snow")
library(snow)
library(tidyverse)
library(stargazer)

#################################################

## In article scenario-based experiment


## Read in data

full_data <- read.csv("scenario_based_data.csv") 


## Inverse scale of difference in reputation so that it is easier to interpret

full_data$diff_rescaled_grep <- -1*full_data$diff_scaled_GREP
full_data$diff_rescaled_drep <- -1*full_data$diff_scaled_DREP

full_data$rescaled_grep <- full_data$Initial_Grep/100

## Summary statistics for Table 2

## separate populations to make summaries easier, will be used for Table 4
modelgood <- subset(full_data, full_data$Good_Rep_Pop==1)
modelmiddling <- subset(full_data, full_data$Middling_Rep_Pop == 1)

mean(modelgood$rescaled_grep)
sd(modelgood$rescaled_grep)

mean(modelmiddling$rescaled_grep)
sd(modelmiddling$rescaled_grep)

## Produce figure 6 

## Use Tidyverse to create a dataframe that includes the mean and standard deviation of 
## difference in expected probability of upholding alliance commitments for each treatment group
## given that they were in the good reputaiton prior population

df <- modelgood%>%
  group_by(Treatment) %>%
  summarize(mean(as.numeric(diff_rescaled_grep)), sd(as.numeric(diff_rescaled_grep))) %>%
  rename(Mean = "mean(as.numeric(diff_rescaled_grep))") %>%
  rename(SD = "sd(as.numeric(diff_rescaled_grep))")


# plot the point plot 


p<-ggplot(df, aes(x=as.factor(Treatment), y=Mean)) + 
  geom_pointrange(aes(ymin=Mean-SD, ymax=Mean+SD), width=.2) 
  
## Add sensible names for the x axis labels
p + scale_x_discrete(name ="Treatment", 
                     labels=c("0" = "No Compliance", "1" = "Partial",
                               "2" = "Full Compliance"))




## Produce figure 7

## Do the same as above but now for the middling population.

df2 <- modelmiddling%>%
  group_by(Treatment) %>%
  summarize(mean(as.numeric(diff_rescaled_grep)), sd(as.numeric(diff_rescaled_grep))) %>%
  rename(Mean = "mean(as.numeric(diff_rescaled_grep))") %>%
  rename(SD = "sd(as.numeric(diff_rescaled_grep))")


# plot the point plot 


p2<-ggplot(df2, aes(x=as.factor(Treatment), y=Mean)) + 
  geom_pointrange(aes(ymin=Mean-SD, ymax=Mean+SD), width=.2) 

p2 + scale_x_discrete(name ="Treatment", 
                     labels=c("0" = "No Compliance", "1" = "Partial",
                              "2" = "Full Compliance"))






## Models for Table 3


model1_global <- lm(full_data$diff_rescaled_grep ~ full_data$Good_Rep_Pop + full_data$Treatment_Middling +
                      full_data$Treatment_Bad +  full_data$Good_Rep_Pop*full_data$Treatment_Middling +
                    full_data$Good_Rep_Pop*full_data$Treatment_Bad+  full_data$age
             + full_data$Female_Dummy + full_data$White_Dummy + full_data$Income + full_data$Party_ID + full_data$Edu)

model1_dyadic <- lm(full_data$diff_rescaled_grep ~ full_data$Good_Rep_Pop + full_data$Treatment_Middling +
                      full_data$Treatment_Bad +  full_data$Good_Rep_Pop*full_data$Treatment_Middling +
                    full_data$Good_Rep_Pop*full_data$Treatment_Bad+  full_data$age
                    + full_data$Female_Dummy + full_data$White_Dummy + full_data$Income + full_data$Party_ID + full_data$Edu)

stargazer(model1_global, model1_dyadic)





## Models for Table 4


modelgood <- subset(full_data, full_data$Good_Rep_Pop==1)
modelmiddling <- subset(full_data, full_data$Middling_Rep_Pop == 1)

modelgood1 <-lm(modelgood$diff_rescaled_grep ~ modelgood$Treatment_Bad + modelgood$Treatment_Middling + modelgood$age
                + modelgood$Female_Dummy + modelgood$White_Dummy + modelgood$Income + modelgood$Party_ID + modelgood$Edu)
modelmiddling1 <-lm(modelmiddling$diff_rescaled_grep ~ modelmiddling$Treatment_Bad + modelmiddling$Treatment_Good + modelmiddling$age
                    + modelmiddling$Female_Dummy + modelmiddling$White_Dummy + modelmiddling$Income + modelmiddling$Party_ID + modelmiddling$Edu)

stargazer(modelgood1,  modelmiddling1)





########################################################################################################################

## Appendices

## Models for Dyadic Reputation Table 5
modelgood2 <- lm(modelgood$diff_rescaled_drep ~ modelgood$Treatment_Bad + modelgood$Treatment_Middling + modelgood$age
                 + modelgood$Female_Dummy + modelgood$White_Dummy + modelgood$Income + modelgood$Party_ID + modelgood$Edu)
modelmiddling2 <- lm(modelmiddling$diff_rescaled_drep ~ modelmiddling$Treatment_Bad + modelmiddling$Treatment_Good + modelmiddling$age
                     + modelmiddling$Female_Dummy + modelmiddling$White_Dummy + modelmiddling$Income + modelmiddling$Party_ID + modelmiddling$Edu)


stargazer(modelgood2, modelmiddling2)


## Attention Check Condition

## Subset the data to just those indivduals that passed at least 4 of 5 attention checks


## Attention Check subset Test 1-- Table 6



attention_check <- subset(full_data, full_data$Attention >= 3) 

model1_attention_check <- lm(attention_check$diff_rescaled_grep~ attention_check$Good_Rep_Pop + attention_check$Treatment_Bad + attention_check$Treatment_Middling + 
                               attention_check$Female_Dummy + attention_check$White_Dummy +
                               attention_check$Edu + attention_check$age + attention_check$Income + attention_check$Party_ID +
                               attention_check$Good_Rep_Pop*attention_check$Treatment_Bad + attention_check$Good_Rep_Pop*attention_check$Treatment_Middling)



model2_attention_check <-lm(attention_check$diff_rescaled_drep~ attention_check$Good_Rep_Pop + attention_check$Treatment_Bad + attention_check$Treatment_Middling + 
                              attention_check$Female_Dummy + attention_check$White_Dummy +
                              attention_check$Edu + attention_check$age + attention_check$Income + attention_check$Party_ID +
                              attention_check$Good_Rep_Pop*attention_check$Treatment_Middling + attention_check$Good_Rep_Pop*attention_check$Treatment_Bad)




stargazer(model1_attention_check, model2_attention_check)


############################################################

## Attention Check subset Test 2 -- Produce Table 7



## Subset data into two treatment populations
modelgood_ac <- subset(attention_check, attention_check$Good_Rep_Pop==1)
modelmiddling_ac <- subset(attention_check, attention_check$Middling_Rep_Pop == 1)

## models with just those that got greater than three on the attention check questions

modelgood1_ac <-lm(modelgood_ac$diff_rescaled_grep ~ modelgood_ac$Treatment_Bad + modelgood_ac$Treatment_Middling + modelgood_ac$age
                   + modelgood_ac$Female_Dummy + modelgood_ac$White_Dummy + modelgood_ac$Income + modelgood_ac$Party_ID + modelgood_ac$Edu)
modelgood2_ac <- lm(modelgood_ac$diff_rescaled_drep ~ modelgood_ac$Treatment_Bad + modelgood_ac$Treatment_Middling + modelgood_ac$age
                    + modelgood_ac$Female_Dummy + modelgood_ac$White_Dummy + modelgood_ac$Income + modelgood_ac$Party_ID + modelgood_ac$Edu)
modelmiddling1_ac <-lm(modelmiddling_ac$diff_rescaled_grep ~ modelmiddling_ac$Treatment_Bad + modelmiddling_ac$Treatment_Good + modelmiddling_ac$age
                       + modelmiddling_ac$Female_Dummy + modelmiddling_ac$White_Dummy + modelmiddling_ac$Income + modelmiddling_ac$Party_ID + modelmiddling_ac$Edu)
modelmiddling2_ac <- lm(modelmiddling_ac$diff_rescaled_drep ~ modelmiddling_ac$Treatment_Bad + modelmiddling_ac$Treatment_Good + modelmiddling_ac$age
                        + modelmiddling_ac$Female_Dummy + modelmiddling_ac$White_Dummy + modelmiddling_ac$Income + modelmiddling_ac$Party_ID + modelmiddling_ac$Edu)

stargazer(modelgood1_ac, modelmiddling1_ac, modelgood2_ac, modelmiddling2_ac)
